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Abstract 

We  present  a  method  for  estimating  the  absolute  frequency  stability  of  N  clocks  separate  from 
a  reference.  The  method  introduced  is  a  modification  of  the  one  proposed  by  Tavella  and  Premoli 
(1993).  After  developing  the  theory  we  apply  the  method  to  atomic  clock  data  gathered  from  the 
USNO. 


INTRODUCTION 

The  estimation  of  absolute  frequency  stability  of  clocks  separate  from  a  reference  clock  is  usually 
required  in  order  to  produce  a  mean  time  scale.  The  usual  method  which  can  be  employed  for 
this  task  is  the  so-called  “N-cornered-hat”  method.  As  first  presented  by  Gray  and  Allan  [1],  this 
method  estimates  the  stability  of  a  fixed  clock  among  an  ensemble  of  N  clocks  by  forming  all 
(N  —  1  )(JV  —  2)/2  triads  of  time  differences  which  include  those  of  the  clock  under  test.  The  reason 
for  triads  is  to  avoid  overdetermination  in  the  estimation  problem.  From  each  of  these  triads  one 
obtains  stability  estimates  for  each  of  the  clocks  in  the  triad  under  the  assumption  that  clocks  are 
uncorrelated  —  these  are  the  “three-cornered-hat”  estimates.  These  estimates  are  then  weighted 
by  a  “triad  uncertainty”  to  yield  a  stability  estimate  for  each  clock.  Another  version,  developed 
by  Barnes  [2],  simultaneously  uses  all  N  —  1  time  differences  with  a  least-squares- type  approach 
to  estimate  the  individual  stabilities.  These  approaches  sometimes  lead  to  stability  estimates  that 
are  negative.  Some  attribute  the  negative  estimate  to  its  uncertainty  [3],  while  others  suggest  the 
possibility  that  the  assumption  of  uncorrelation  is  not  always  valid  [4] ;  both  are  probably  true  to 
a  significant  extent. 

Recently,  Tavella  and  Premoli  [4]  developed  an  approach  which  avoids  the  problem  of  negative 
stability  estimates  by  allowing  the  possibility  of  correlations  among  the  ensemble  of  clocks.  Their 
approach  is  consistent  and  formally  equivalent  to  the  three-cornered-hat  method  when  clocks  are 
uncorrelated,  that  is,  both  methods  produce  the  same  result  in  this  situation.  Their  method, 


69 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

DEC  1998 


2.  REPORT  TYPE 


3.  DATES  COVERED 

00-00-1998  to  00-00-1998 


4.  TITLE  AND  SUBTITLE 

Estimating  Frequency  Stability  and  Cross-Correlations 


6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

U.S.  Naval  Observatory ,3450  Massachusetts  Ave 
NW, Washington, DC, 20392 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROIECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 


11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

See  also  ADA415578.  30th  Annual  Precise  Time  and  Time  Interval  (PTTI)  Systems  and  Applications 
Meeting,  Reston,  VA,  1-3  Dec  1998 


14.  ABSTRACT 

see  report 


15.  SUBIECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 


a.  REPORT 

unclassified 


b.  ABSTRACT 

unclassified 


c.  THIS  PAGE 

unclassified 


17.  LIMITATION  OF 

18.  NUMBER 

ABSTRACT 

OF  PAGES 

Same  as 

13 

Report  (SAR) 

19a.  NAME  OF 
RESPONSIBLE  PERSON 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


however,  will  not  produce  a  negative  variance  estimate  even  when  a  three-cornered-hat  approach 
would.  Since  the  problem  of  estimating  absolute  frequency  stability  is  underdetermined,  they 
impose  the  condition  that  if  correlations  exist  they  should  be  “small.”  They  propose  a  condition 
which  amounts  to  finding  the  smallest  correlations  possible  which  keep  a  certain  matrix  positive- 
definite.  This  will  be  explained  in  greater  detail  below.  Since  Tavella  and  Premoli’s  work  [4] 
was  done  in  the  framework  of  three  clocks,  a  weighted  triad  uncertainty  approach  can  be  used  to 
improve  the  estimates  if  there  are  more  than  three  clocks  in  the  ensemble. 

In  a  later  paper  [5]  Tavella  and  Premoli  refined  their  analysis.  They  show  that  when  the  number 
of  clocks  in  a  comparison  increases,  the  amount  of  arbitariness  in  determining  absolute  frequency 
stability  decreases. 

The  goal  of  this  paper  is  to  generalize  the  methods  of  [4]  and  [5]  and  to  estimate  absolute  frequency 
stability  for  N  clocks  (N  >  3).  The  results  of  this  analysis  will  be  discussed  and  applied  to  atomic 
clock  time  difference  data  from  the  U.  S.  Naval  Observatory  (USNO). 


NOTATIONS  AND  DEFINITIONS 

Let  xl  —  {x\  :  fc  >  1}  for  i  —  1,  •  ■  ■ ,  JV  —  1  be  the  time  differences  of  clock  i  with  respect  to  a 
fixed  reference  clock,  say  clock  TV,  which  is  sampled  every  to  seconds.  We  let  XJ  —  {X3k  :  k  >  1} 
for  j  =  1  ,•••  ,N  represent  the  time  differences  with  respect  to  a  noiseless  “ideal”  clock,  so  that 
xl  =  X1  —  XN. 


Define  the  fractional  frequencies: 


xL. ,  -  xl 


S  k  >  1 


and  the  process  of  averaging  such  fractional  frequencies  of  order  m  >  1: 

-i  f  \  y(k—l)m+l  V(fc-l)m+2  y'km  1  j 

yU™ro)  =  — — —  =  -  E  y\k-i)m+i •  (2) 

Notice  that  (2)  reduces  to  (1)  for  m  —  1.  We  define  Yk  and  Yk(rriTo)  by  replacing  x  with  X 

appropriately.  As  a  measure  of  time  domain  frequency  stability  we  use  the  Allan  variance.  The 
Allan  variance  of  clock  i  referenced  to  clock  N  is  defined  as 

sn(r)  =  7}((y2(r)  -  yl(T))2)  (3) 

where  {)  denotes  mathematical  expectation.  We  assume  here,  of  course,  that  the  process  of  averaged 

fractional  frequencies  is  stationary  and  ergodic  so  that  the  definition  (3)  is  well-defined.  Most 
authors  use  the  notation  of  N(r)  to  represent  the  quantity  in  (3),  but  we  will  use  the  notation  set 
forth  in  the  works  [4]  and  [5].  For  what  follows  we  need  to  define  the  Allan  covariance  of  clocks  i 
and  j  referenced  to  clock  N  to  be 


i(T)  =  o<(i/2(t)  *  y\{T))(yi(r)  -  yi(r))). 
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This  statistic  appears  in  [6]  with  different  notation.  It  is  worth  mentioning  at  this  point  that  a 

covariance  reduces  to  a  variance  for  i  —  j;  therefore  we  will  refer  to  variances  as  covariances  when  no 
confusion  can  arise.  Furthermore,  it  is  clear  from  the  definition  (4)  that  Sij(r)  —  Sji(r).  However, 
all  these  quantities  are  usually  never  available  in  practice  and  can  only  be  estimated  through  a 
large  stretch  of  data  (at  least  for  r  =  mr0  and  integer  m  >  1)  by  the  samples 


1 


71— 2m 


— 


£  (yi+AmTo)  -  yHmro))(yi+i(mT°'>  ~y£(mr &)) 


v  2  (n  —  2m)  ^ 

where  n  is  the  data  length.  The  goal  of  this  paper  is  estimate  the  dereferenced  quantities: 

1 


(6) 


and 


1 


Tij  =  2<(n*(r)  -  Y{{r)){Yi(r)  -  Y(( r))).  (7) 

Again,  some  authors  use  the  notation  a f(r)  to  represent  the  quantity  shown  in  (6).  The  associated 
sample  quantity  is  defined  similar  to  (4): 


1 


7i—2  m 


ra 


n  —  2  m 


53  (Yk+  l(mr0)  -  Yk(mTo))(Yk+l(mro)  -  Yk(mro)) 


k= 1 


By  noting  that  y’k  =  Y£  -  Yj^  we  deduce  from  (4)- (8)  the  following  relations; 

Sij  =  rtj  +  r^N  —  r  iN  —  rjN 

and 

Sij  —  Tij  +  rjvjv  “  fiN  —  ?jN- 
If  IV  —  3  and  all  fij  —  0  for  i  ^  j  then  the  usual  three-cornered-hat  estimates  fall  out  of  (10): 

fll  =  Sll  -  Sl2 
?22  ~  S22  ~  S 12  ■ 

^33  —  5 12 


(8) 

(9) 

(10) 

(11) 


This  was  noted  in  references  [4]  and  [6]. 

THE  TAVELLA-PREMOLI  APPROACH 

Let’s  first  outline  the  approach  followed  by  Tavella  and  Premoli  [4].  Suppose  we  are  interested  in 
obtaining  estimates  of  absolute  frequency  stabilities  for  three  clocks,  say  clocks  1,  2,  and  3.  That 
is,  we  wish  to  estimate  the  matrix 


R  = 


r  12 

nz 

7*12 

r22 

r23 

7*13 

^23 

rZZ 

(12) 
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Notice  the  matrix  R  is  inherently  symmetric.  Now  consider  the  matrix 


S  = 


■511  S\2 

$12  $ 22 


S  is  symmetric  and  we  know  a  priori  S  is  also  positive-definite.  The  relation  (9)  can  be  rewritten 
to  give  the  relations: 

ni  =  s  ii  —  r33  +  2  ri3 

T22  —  $ 22  —  ^33  +  2f23 

T\2  ~  s12  —  r33  +  rl3  +  r23. 


This  means  that  all  the  entries  of  R  can  be  written  as  functions  of  the  elements  r  —  (ri3,  r23)  and 
the  reference  stability  r33.  To  find  estimates  for  the  entries  of  R ,  Tavella  and  Premoli  proposed  the 
following.  We  should  find  values  of  covariances  ri2,ri3  and  r23,  which  are  small  in  some  sense,  but 
keep  the  matrix  R  positive-definite.  To  understand  in  what  sense  they  mean  small  define 

G(ri2,ri3,r23)  =  r%2  +  r?3  +  rf3.  (13) 


Clearly,  the  minimum  value  of  G  is  0  exactly  when  ri3  =  r23  =  7*12  =  0,  and,  given  (9),  this  implies 
7-33  =  sl2.  That  is,  we  arrive  at  the  three-cornered-hat  estimates  (11): 


^3 


-corner  ~ 


—  s12  0  0 

0  S22  —  S12  0 

0  0  5i2 


(14) 


As  long  as 


«n  >  $12 

$22  >  sn  (15) 

a  12  >  0 


we  have  valid  estimates  of  stability.  However,  if  one  of  the  diagonal  elements  in  (14)  is  negative,  we 
obtain  a  negative  stability  estimate.  This  can  be  ruled  out  if  we  insist  that  the  estimated  matrix 
R  be  positive-definite.  Tavella  and  Premoli  are  able  to  show  that  the  function 

1  (p 

H (ri3,  r23,  r33)  —  r33  —  (ri3  —  r33,  r23  —  r33)5~  (ri3  —  r33,  r23  -  r33)  (16) 

is  positive  if  and  only  if  R  is  positive-definite,  T  denotes  transpose  here.  Therefore  they  suggest 
the  rij  should  be  chosen  in  a  way  that  minimizes  the  following  expression: 

*•12  +  rh  +  r23  /l7x 

H{ris,r23,r33) 

where  the  minimum  is  taken  over  all  values  ri3, r23, r33  for  which  H  is  positive  (i.e.,  R  is  positive- 
definite).  Notice  that  if  H(0, 0,  .S'12)  >  0  then  the  minimum  of  (17)  is  0  and  we  again  achieve 
the  three-cornered-hat  estimates.  This  procedure,  however,  will  not  lead  to  a  negative  stability 


estimate.  To  see  this  just  notice  that  the  minimization  occurs  over  a  convex  region,  namely  the 
elliptic  paraboloid  region 

P  =  {{x,y,z)  ■  H{x:y,z)  >  0}. 

The  boundary  of  P,  6P ,  contains  all  the  points  where  the  function  H  identically  equals  0.  Since  the 
objective  function  (17)  is  convex,  it  is  well  known  that  the  minimizer  of  this  problem  is  unique.  The 
minimizer  (r*J3,  r^,  t~33)  of  (17)  will  always  satisfy  H (r^,  r^,  r^)  >  0  since  points  on  the  boundary 
of  P  yield  undefined  values  in  (17).  Now  H  >  0  implies  R  is  positive-definite.  The  result  follows 
from  the  fact  that  positive-definite  matrices  have  positive  diagonal  elements.  The  minimization  of 
(17)  was  carried  out  explicitly  for  IV  =  3  clocks  [4]  and  its  solution  was  formulated  through  the 
zeroes  of  a  sixth-degree  polynomial. 

The  implementation  of  this  method  does  show  some  subtle  inadequacies  however.  As  noted  in  the 
introduction,  a  three-cornered-hat  approach  may  give  a  negative  stability  estimate.  Although  the 
Tavella-Premoli  approach  will  not  produce  such  a  negative  value,  it  may  produce  an  optimistically 
small  estimate  of  stability.  Let’s  see  why  this  may  happen.  It  is  easy  to  show  that  if  the  conditions 
(15)  are  all  satisfied  then  (0, 0,  Su)  €  P.  Therefore,  if  ,§12  is  arbitrarily  small  and  positive,  our 
estimate  for  r 33  will  also  be  arbitrarily  small.  Similarly,  if  sn  <  S22  and  S12  is  close  to  sqi ,  our 
estimate  for  rq  may  also  be  optimistically  small.  These  effects  have  appeared  in  time  difference 
data  gathered  at  the  USNO. 

This  optimism  can  be  attributed  to  the  rather  large  domain  of  admissible  values  in  (17).  The 
admissible  region  is  an  elliptic  paraboloid  in  Tf3 .  If  more  clocks  are  in  our  comparison  then,  as 
shown  in  [5],  this  admissible  region  is  substantially  reduced  and  points  that  would  be  admissible 
when  considered  through  triads  would  be  disallowed  when  considered  in  a  multiple  comparison.  We 
can  use  this  substantial  reduction  in  admissible  values  to  generalize  the  Tavella-Premoli  scheme. 
We  will  take  this  up  next. 


A  MODIFIED  APPROACH 


Suppose  we  have  time  differences  from  N  -  1  clocks  with  a  fixed  reference  clock  TV  (TV  >  3)  and  we 
wish  to  estimate  the  stabilities.  We  could,  of  course,  apply  the  method  of  Tavella  and  Premoli  [4]  to 
triads  of  clocks  and  weight  appropriately.  As  noted  earlier,  this  approach  may  produce  optimistic 
values  of  stability. 


Let’s  consider  the  following  approach.  In  analogy  to  the  method  proposed  by  [4],  we  suggest  an 
obvious  modification  to  (17)  and  find  the  values  of  covariances  that  minimize  the  following  Tavella- 
Premoli  function: 


\  ,r?. 
■n<j  ij 


H2{r\N,  -  •  ■  ,Tnn) 


(18) 


where 


■  •  ■  ,tnn)  —  r^N  ~  (n n  -  Dvtv,  •  •  •  -  rNN)S  1(r1/v  -  rjvjv,  ■  •  •  ,ruv  -  rNN)T 

the  superscript  T  denotes  transpose  and  the  minimization  is  over  those  points  r,Ar  which  keep  H 
positive.  Here,  the  function  H  >  0  if  and  only  if  the  matrix  R  is  positive-definite  [5].  Dividing 
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by  H  above  as  mentioned  earlier  keeps  the  points  away  from  the  boundary  of  the  admissible 
region  and,  consequently,  will  keep  the  matrix  R  positive-definite.  Notice  the  function  (18)  is 
now  a  function  of  the  N  variables  (r^v,  •  •  • ,  r\N)  as  well  as,  of  course,  the  .s,j.  Also  notice  that 
function  (18)  has  H  squared  in  the  denominator.  This  squaring  is  suggested  in  order  to  keep  the 
minimization  problem  scale-invariant.  This  will  also  be  helpful  in  order  to  do  the  minimization  from 
a  numerical  point  of  view.  The  problem  of  minimizing  (18)  is  a  constrained  minimization  problem 
since  we  are  only  allowing  values  of  r y  which  keep  the  matrix  R  positive-definite,  that  is,  those  rt] 
for  which  H(r\Ni  ■  ■  • ,  rNN)  >  0.  Prom  the  same  convexity  considerations  the  minimization  problem 
(18)  has  a  unique  solution.  Usually  constrained  minimization  problems  are  difficult  to  solve,  but 
one  can  apply  numerical  techniques  if  care  is  taken  in  scaling  and  choice  of  initial  data.  We  used  a 
conjugate  gradient  method  to  produce  the  minimizer.  As  far  as  the  choice  of  initial  data  we  chose 

riN  —  0  for  i  <  N, 

1  1 

rNN  =  —  — 

s*  =  (i,---,i)s-Hi,---,i)T  >  0  from  the  positive-definiteness  of  the  matrix  S  and  thus  S  1 . 
The  factor  \  above  is  used  to  force  the  initial  data  to  lie  within  the  constraint  (using  convexity  of 
admissible  values  here).  This  choice  of  initial  data  conforms  with  our  belief  that  clocks  are  close 
to  uncorrelated.  Of  course,  from  the  uniqueness  of  the  solution  any  reasonable  point  we  choose 
initially  will  converge  to  the  minimizer  and  this,  in  fact,  is  observed.  We  should  mention  that 
problems  can  arise  if  the  initial  datum  is  too  close  to  the  boundary  of  the  admissible  region.  If  this 
is  the  case  a  conjugate  gradient  direction  may  lead  you  outside  the  admissible  region. 

It  is  interesting  to  note  that  if  clocks  are  uncorrelated  (rtJ  —  0  for  all  i  ^  j)  the  minimization  of 
(18)  leads  to  the  estimate 

—  ^  j  $ij 

rNN  ~  (N  —  1){N  —  2)  ’ 

that  is,  the  straight  average  of  the  values  that  one  would  obtain  by  performing  the  usual  Tavella- 
Premoli  procedure  on  all  triads  of  clocks. 

One  implicit  point  we  have  stressed  in  this  paper  is  the  possibility  of  the  existence  of  clock  noise 
correlations.  In  order  to  convince  the  reader  that  clock  correlations  are  a  reality  we  have  developed 
a  preliminary  statistical  test  which,  although  heuristic,  may  be  used  to  determine  if  clock  noise 
correlations  exist  between  clocks  in  the  comparison. 


A  STATISTICAL  TEST  OF  CORRELATION 


Tavella  and  Premoli  have  shown  [5]  that  if  rl3  —  0  for  all  i  7^  j  then  the  matrix  S  of  “true”  Allan 
covariances  has  the  form 


S  = 


r  11  -  tnn 
uvjv 
rNN 


rNN  tnn 

T22  -  rNN  rNN 

rNN  r33  —  r  nn 


rNN 

rNN 

rNN 


(19) 


rNN  rNN  rNN 


rN-i,N-i  -  rNN 
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As  mentioned  in  the  previous  section,  the  matrix  5  can  only  be  estimated  by  the  matrix  S  =  (sij). 
Thus,  if  the  clocks  axe  uncorrelated  then  the  stJ  for  i  7^  j  are  all  estimating  the  reference  stability 
r/vjv.  If  the  data  length  n  is  large,  one  should  expect  f{j  «  0  for  i  ^  j.  If  this  is  the  case  then 
from  (10)  and  (19)  the  statistic  Sij  «  frviv  for  i  7^  j.  Now  since  the  distribution  of  r  jvjv  can 
be  shown  to  behave  as  a  chi-square  with  some  number  of  degrees  of  freedom  (we  consider  here 
overlapping  sample  estimates),  which  in  general  varies  with  the  noise  type  present  [3],  we  can  test 
if  the  off-diagonal  terms  of  S  actually  differ  up  to  some  statistical  significance.  The  problem  can 
be  formulated  as  follows: 

Consider  as  the  null  hypothesis  all  clocks  are  uncorrelated.  We  would  like  to  test  whether  the 
alternative  hypothesis  that  some  pair  of  clocks  exhibit  correlations  is  true.  So  we  would  like  to 
test  for  the  simultaneous  equality  of  the  off-diagonal  elements  of  the  matrix  S.  This  is  a  multiple 
comparison  problem  whose  analysis  can  be  quite  difficult.  We  instead  will  be  content  to  work  with 
pair-wise  comparisons  of  the  off-diagonal  elements  Sjj .  For  this  problem  we  test  whether  the  ratio 


for  different  pairs  of  clocks  i  7^  j  and  1’  yf  j'  is  significantly  different  from  1.  Now  since  the  chi- 
squares  all  have  the  same  degrees  of  freedom,  d,  for  a  fixed  integration  time  r,  the  test  statistic 
/  has  the  i^-distribution,  that  is,  the  ratio  of  two  chi-squares  with  d  degrees  of  freedom  each. 
Assume  that  for  fixed  r  >  0  all  the  off-diagonal  terms  of  S  are  positive.  For  90%  confidence  we 
choose  to  reject  the  null  in  favor  of  the  alternative  if  either 

/  >  Fd,d,.9 5  or  /  <  Fd,d,. 05 

where  Fd,d,a  represents  the  ath  percentile  of  the  F  distribution.  Since  all  estimates  §ij  share  the 
same  degrees  of  freedom,  the  problem  is  simplified  considerably  since  we  only  need  to  check  if  the 
statistic 

j-* 

f  =  max  ~- 

Si’  j’ 

satisfies  /*  >  Fd,d,, 95.  The  above  argument  is  only  heuristic  and  the  resulting  statistical  test  is  not 
entirely  substantiated  since  significant  accidental  covariances  for  increasingly  smaller  sample  sizes 
can  corrupt  the  distributional  properties  of  the  /-statistic,  that  is,  if  n  —  2 m  is  not  “large”  then  the 
values  fki  for  k  7^  l  in  (10)  may  not  be  close  to  zero  and  can  bias  the  -statistic  away  from  f^N- 
However,  at  least  for  very  large  sample  sizes  n  and  relatively  small  m,  the  above  approximation 
seems  reasonable.  Of  course,  a  better  statistical  test  can  be  achieved  if  one  can  characterize  the 
distribution  of  the  Allan  covariance  statistics. 

We  applied  the  above  analysis  to  a  group  of  four  cesium-beam  frequency  standards  rereferenced 
to  a  fifth  cesium-beam  standard  (see  Table  1).  The  time  differences  had  data  length  n  —  167,513. 
We  made  the  assumption  that  for  the  integration  times  observed  the  dominant  noise  type  was 
white  frequency  modulation  and  computed  the  appropriate  degrees  of  freedom  for  the  overlapping 
estimates: 

r3(n  —  1)  2(n  — 2)n  4m2 

L  2m  n  J  4m2  -H  5 ' 

We  found  the  corresponding  F  quantile  and  applied  the  above  arguments.  We  noticed  that  although 
the  off-diagonal  entries  of  S  may  seem  to  be  the  same,  at  least  for  short  integration  times  (from 


20  seconds  out  to  about  320  seconds)  we  were  able  to  reject  the  null  hypothesis  and  conclude  that 
statistically  significant  correlations  exist.  This  does  not  say  that  a  correlation  is  ‘large”;  indeed, 
rij  can  be  relatively  quite  small  (by  orders  of  magnitude)  when  compared  to  the  diagonal  entries  ra 
and  fjj .  However,  some  small  covariance  seems  quite  reasonable  when  we  consider  that  some  of  the 
clocks  in  the  above  analysis  were  in  the  same  environmental  chamber.  An  exhaustive  treatment  of 
the  above  ideas  would  certainly  be  interesting  from  both  a  practical  and  theoretical  point  of  view. 


5(20)  - 


5(40)  = 


5(320)  = 


7.10826#  -  24 
3.81328 E  -  24 
3.79768#  -  24 
3.79259#  -  24 

3.22437#  -  24 
1.69300# -24 
1.69913#  -  24 
1.69097#  -  24 

3.55895#  -  25 
1.82808#  -  25 
1.85889#  -  25 
1.84763#  -  25 


3.81328#  -  24 
7.95851#  -  24 
3.82652#  -  24 
3.83888#  -  24 

1.69300#  -  24 
3.42756#  -  24 
1.69388#  -  24 
1.70435#  -  24 

1.82808#  -  25 
3.65834#  -  25 
1.89250#  -  25 
1.85393#  -  25 


3.79768#  -  24 
3.82652#  -  24 
7.89671#  -  24 
3.82095#  -  24 

1.69913#  -  24 
1.69388#  -  24 
3.58596#  -  24 
1.71195# -24 

1.85889#  -  25 
1.89250# -25 
4.04481#  -  25 
1.89242#  -  25 


3.79259#  -  24  ' 
3.83888#  -  24 
3.82095#  -  24 
6.99711#  -  24 

1.69097#  -  24 
1.70435#  -  24 
1.71195#  -24 
3.15194#  -  24 

1.84763#  -  25 
1.85393#  -  25 
1.89242#  -  25 
3.55988#  -  25 


Table  1. 


Averaging  Time 

F-statistic 

/*-statistic 

20  sec. 

1.009893 

1.012205 

f*  >  F:  can  reject  in  favor  of 
alternative. 

40  sec. 

1.010690 

1.012407 

/*  >  #:  can  reject  in  favor  of 
alternative. 

320  sec. 

1.026667 

1.035239 

/*  >  #:  can  reject  in  favor  of 
alternative. 

The  S  matrices  above  were  computed  for  r  =20,  40  and  320  seconds.  Table  1  shows  the  results  of 

the  statistical  test  for  these  integration  times. 


RESULTS 

We  have  applied  the  technique  described  in  the  previous  section  to  time  difference  data  obtained 
from  atomic  clocks  at  the  U.S.  Naval  Observatory  (USNO).  USNO  has  a  large  ensemble  of  clocks 
consisting  of  both  commercial  cesium-beam  standards  and  active  hydrogen  masers.  Data  from  both 
types  of  clocks  were  considered  separately. 


Figures  la  and  lb  show  sigma-tau  plots  generated  by  a  three-cornered-hat  and  Tavella-Premoli 
analysis  respectively.  For  comparison,  both  analyses  were  carried  out  on  cesium-beam  standards 
labelled  229,  260,  and  253  over  a  period  of  five  years.  The  three  time  series  were  rereferenced  to 
one  of  the  clocks  (253)  to  avoid  the  obvious  problem  of  correlation  between  the  time  scales.  We 
used  overlapping  estimates  of  Allan  covariance  over  all  r,  instead  of  just  r  =  2mro  (integer  m  >  1) 
where  r0  is  the  sampling  time. 

There  are  several  things  to  notice  about  these  plots.  First,  both  approaches  coincide  where  the 
three-cornered-hat  approach  produces  positive  stability  estimates,  so  that  the  primary  difference 
between  the  two  is  that  the  Tavella-Premoli  scheme  is  able  to  produce  an  ”  extension”  of  the 
reference  stability.  Second,  by  plotting  the  data  for  all  r,  a  ”  pathological”  behavior  is  apparent  in 
both  approaches.  There  are  two  ’’dips'1  in  each  graph:  one  at  60  days  and  the  other  at  about  590 
days.  Neither  of  these  dips  are  likely  to  be  physical, since  cesium-beam  standards  integrate  down 
with  a  1/yff  dependence  until  they  reach  their  flicker  floor  or  begin  to  drift  significantly.  One 
approach  to  this  situation  might  be  to  add  another  clock  to  the  ensemble  and  use  the  additional 
three  triads.  However,  even  if  a  straight  average  of  the  resulting  estimates  coming  from  each  triad 
were  used,  one  would  still  find  that  the  dip  creates  a  large  bias  in  the  result. 

Instead,  the  extra  clock  can  be  used  with  the  modified  approach  outlined  in  the  previous  section. 
The  corresponding  frequency  stability  plot  is  shown  in  Figure  2.  This  approach  is  in  close  agreement 
with  the  other  two  until  r  «  30  days.  At  this  point,  the  first  spike  apparent  in  the  other  approaches 
is  smoothed  out.  Most  importantly,  the  new'  approach  now  maintains  the  expected  1  jyfr  behavior 
at  this  point  in  the  data.  The  second  spike  is  smoothed  out  as  well. 

We  performed  this  modified  approach  on  the  data  with  two  more  clocks  added  in.  The  results  are 
shown  in  Figures  3a  and  3b.  Again  the  absence  of  non-physical  dips  is  apparent  and  agreement 
between  this  and  results  shown  in  Figure  2  is  good.  In  all  plots  we  noticed  that  one  of  the  cesium 
clocks  appeared  to  be  much  worse  than  the  others.  Later  we  were  able  to  show  that  this  particular 
clock  had  a  significant  drift  component  that  wasn’t  compensated  for. 

We  also  applied  the  modified  approach  to  active  hydrogen  maser  data.  We,  unfortunately,  were 
not  able  to  extract  such  a  long  stretch  of  data  as  with  the  cesiums.  The  resulting  plot  is  in  Figure 
4. 

We  would  also  like  to  mention  that  the  modified  approach  did  not  differ  significantly  with  the 
manufacturer  specifications,  and,  in  most  cases,  conformed  to  them  closely.  Although  the  results 
of  this  analysis  are  by  no  means  complete,  we  believe  this  method  shows  significant  promise  in 
achieving  estimates  of  absolute  frequency  stability. 


CONCLUSIONS 

The  modified  approach  for  estimating  frequency  stability  generally  performs  well:  it  gives  essentially 
the  same  results  for  short-term  stabilities  (r  <  30  days)  and  non-pathological  behavior  for  long¬ 
term  stabilities  (r  ~  1  year).  We  are  encouraged  by  these  preliminary  results  and  are  hopeful  that 
the  method  can  be  used  to  give  an  accurate  representation  of  long-term  stability  in  atomic  clocks. 
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Fig.  3:  Plot  of  absolute  frequency  stability  using  modified  approach  with  six  clocks  -  stabilities 
for  all  six  clocks  shown  (cf.  Fig.  2).  The  analysis  which  gave  this  plot  was  performed  on  dyadic 
averaging  times  only  as  opposed  to  all  averaging  times  shown  in  Fig.  2. 
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Fig.  4:  Plot  of  absolute  frequency  stability  using  modified  approach  with  four  active 
masers  sampled  at  20-second  intervals.  Analysis  performed  with  dyadic  averaging 


hydrogen 

times. 


Questions  and  Answers 


PATRIZIA  TAVELLA  (IEN):  I  am  happy  that  you  are  continuing  with  this  work  and  you  took  the  time  to 
read  my  papers.  Thank  you  for  the  quotation. 

I  think  that  maybe  the  optimizer  function  can  he  different  in  the  short  and  long  term.  Maybe  the  long 
correlation  is  reasonable  in  the  short  term,  but  not  in  the  long  term.  Maybe  we  could  investigate  if  the 
possible  function  to  be  minimized  should  be  different  in  long  and  short  term. 

I  would  also  like  to  investigate  the  effect  of  the  uncertainty'  of  the  variance  evaluation  on  the  positive 
definitiveness  of  the  matrix,  because  since  the  evaluation,  I  am  uncertain  due  to  the  limited  number  of 
measurements.  Maybe  the  region  can  change. 

FRED  TORCASO  (USNO):  I  think  moving  in  that  direction  would  be  trying  to  get  a  better  understanding 
of  the  distribution  of  tire  sample  Allan  variance  for  exactly  those  integration  times.  That  is  actually  a 
difficult  problem.  There  have  been  a  lot  of  papers  appearing  which  try  to  estimate  the  distribution  of  the 
sample  Allan  variance  for  integration  times  of,  I  think,  eight  samples.  It  is  something  ridiculously'  hard. 

It  looks  like  there  is  some  evidence  for  large  integration  times  in  many  sample  estimates  that  a  gaussian 
approximation  is  probably  a  good  one.  Maybe  this  would  be  a  possible  way  to  go. 

MARC  WEISS  (NIST):  Have  you  thought  about  how  to  compute  uncertainty  for  Allan  variance  when  you 
have  the  n-comered-hat  technique?  In  other  words,  there  is  an  uncertainty  for  the  Allan  variance  when 
you  look  at  the  distribution  just  because  you  have  a  sample  variance. 

FRED  TORCASO:  So  you  like  to  put  error  bars  around  the  absolute  values? 

MARC  WEISS:  In  the  presence  of  an  n-comered  hat  where  you  are  estimating  a  variance  through  the 
other  clocks,  I  think  there  is  an  additional  uncertainty  because  of  that.  I  think  the  relative  stability  of  the 
clocks  will  come  in  as  well.  In  other  words,  if  you  look  at  one  very  good  clock  and  you  have  ten  clocks 
that  are  ten  times  worse,  you  can  not  see  the  good  clock  -  at  least  in  theory  you  should  not  be  able  to. 

I  do  not  know'  how  to  do  it.  I  am  wondering  if  you  have  looked  at  that? 

FRED  TORCASO:  I  actually  have  not  looked  at  ways  of  obtaining  error  bars  on  estimates  of  absolute 
stability.  I  am  looking  into  the  possibility  of  comparing  very  quiet  clocks,  for  instance,  hydrogen  masers, 
over  a  short  time  scale  compared  with  the  high-performance  cesium  beam  clocks  that  we  have  at  USNO. 

A  similar  approach  to  the  one  I  described  may  work  if  we  rescale  this  co-variance  matrix  that  I 
introduced,  “R”,  to  a  correlation  matrix.  Then  all  the  off-diagonal  terms  of  the  matrix  are  of  the  same 
order.  That  will  help  in  the  analysis,  but  I  have  not  done  that. 


